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1 . I nt roduct i on 


In  geodesy  and  geophysics  we  frequently  meet  with  the 
situation  that  a model  defined  by  a set  of,  say,  N parameters 
is  to  be  determined  from  a smaller  number  n < N of  observations. 

As  an  example,  the  internal  structure  of  the  earth  may  be 
defined  by  a set  of  N parameters  describing  the  density,  the 
rigidity,  and  the  compressibility  of  the  earth  as  a function  of 
depth.  The  n observations  comprise  velocities  of  seismic  sur- 
face waves,  together  with  the  mass  and  the  polar  moment  of  inertia 
of  the  earth.  If  the  model  for  the  earth's  internal  structure  is 
to  be  realistic,  then  N will  be  large  and  n < N . 

We  thus  have  n < N equations  for  N unknowns,  which  is 
obviously  an  underdetermined  problem  admitting  an  infinite  number 
of  possible  solutions.  Using  standard  mathematical  terminology 
(Lanczos,  1964),  we  have  an  improperly  posed  problem.  (A  problem 
is  properly  posed  if  it  has  a unique  solution  that  depends  con- 
tinuously on  the  data.) 

Originally,  the  equations  expressing  the  data  xi  as 
functions  of  the  model  parameters  sr  will,  in  general,  be  non- 
1 i near : 


;i  ■ Vvs2 V 


i - 1 , 2 n . 


(1-1) 


By  a suitable  application  of  Taylor's  theorem  it  is  usually 
possible  to  approximate  these  equations  by  linear  ones: 


N 

: . = }'  a . s 

i ir  r 


(1-2) 


or  in  matrix  notation: 


2 


The  formal  solution  of  this  system  of  linear  equations 
may  be  written  as 

£ = A_1x  . (1-4) 

If  A were  a regular  square  matrix,  then  A-1  would  be  the  or- 
dinary inverse  matrix  of  A . In  our  underdetermi ned  case,  how- 
ever, A-1  must  be  understood  in  the  sense  of  general ized  matrix 
inverses  (cf.  Bjerhammar,  1973;  Rao  and  Mitra,  1971). 

At  any  rate,  the  solution  of  (1-1)  or  (1-3)  may  be  con- 
sidered as  an  inversion  of  these  equations  with  respect  to  the 
parameters  s r , which  accounts  for  the  name,  geophysical  inverse 
problems. 

Another  typical  example  of  an  "improperly  posed"  inverse 
proolem  is  the  determination  of  subsurface  mass  distributions 
which  produce  a given  anomalous  gravity  field  at  the  earth's  sur- 
face. This  problem  is  sometimes  called  an  inverse  problem  of  po- 
tential theory  (Lavrentiev,  1967;  Burkhard  and  Jackson,  1976). 

The  determination  of  the  earth's  external  gravitational 
field  from  geodetic,  gravimetric  and  satellite  data  may  also  be 
considered  as  an  inverse  problem  that  is  mathematically  quite 
similar  to  the  determination  of  the  internal  structure  of  the 
earth  from  seismic  and  other  data. 

This  geodetic  inverse  problem  is  likewise  underdetermined. 
The  external  gravitational  field  requires  for  a complete  de- 
scription an  infinite  number  of  parameters,  for  instance,  the  set 
of  all  coefficients  in  the  expansion  of  the  external  gravitational 
potential  in  spherical  harmonics.  This  infinite  number,  N = ■»  , 
of  parameters  is  to  be  determined  from  a finite  number  n of 
observati ons . 

Even  in  the  seismic  inverse  problem  it  is,  at  least  theo- 
retically, appropriate  to  take  N = °°  if  we  wish  to  admit  rea- 


sonably  general  functions  for  density,  rigidity,  and  compressi- 
bility because  it  cannot  be  assumed  a priori  that  such  functions 
depend  on  a finite  number  of  parameters  only. 

Thus,  in  general,  the  space  of  parameters  will  be  infinite- 
dimensional  rather  than  N -d i mens i ona 1 . In  other  words,  the  proper 
general  setting  for  (linear)  geodetic  and  geophysical  inverse- 
problems  will  be  i nf i n i te-d i men s i ona I Hilbert  space.  This  was 
pointed  out  by  Krarup  (1969)  for  the  geodetic  case  and  by  Backus 
(1970)  for  geophysical  inverse  problems. 

The  geodetic  inverse  problem,  the  determination  of  the 
external  gravitational  field  from  data  of  different  kind,  is  usu- 
ally solved  by  least-squares  collocation.  This  technique  has  many 
features  in  common  with  other  geophysical  inversion  methods.  It 
may,  therefore,  be  of  interest  to  compare  these  techniques  and  to 
exhibit  some  cross-connections. 

We  shall  also  discuss  least-squares  collocation  from  the 
point  of  view  of  analytically  representing  the  external  gravi- 
tational field  by  a linear  combination  of  suitable  simpler  har- 
mon i c funct  i ons  . 

The  subject  of  the  present  report  is  purely  conceptual, 
aiming  at  a better  understa nd i ng  of  least-squares  collocation  by 
considering  it  in  is  relation  to  other  methods,  no  new  compu- 
tational formulas  will  be  derived.  Still,  this  paper  might  be 
useful  as  a contribution  to  the  present  discussion  on  the  con- 
ceptual foundations  of  least-squares  collocation. 


2 . Systems  of  Linear  Equations 

Let  us  assume  that  the  geophysical  or  geodetic  inversion 
problem  has  already  been  linearized,  so  that  it  reduces  to  the 
solution  of  a system  of  linear  equations  of  form  (1-3), 
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£s  = £ , ( 2 - 1 ) 

£ being  the  vector  of  N parameters  to  be  determined,  £ de- 
noting the  vector  of  n observations,  arid  13  being  a given 
nxN  matrix  of  coefficients. 

Assuming  n < N , we  have  an  u nd erdetermi ned  problem.  The 
number  N of  parameters  may  be  finite  or  infinite.  For  N = «= 
we  must  presuppose  that  the  occurring  sums  from  1 to  N , which 
are  now  infinite  series, converge;  otherwise  the  formal  opera- 
tions are  the  same  as  for  a finite  N . 

Equations  of  type  (2-1)  may  be  formulated  for  seismic 
inversion  problems  (Knopoff  and  Jackson,  1973),  for  gravity  inter- 
pretation problems  (Burkhard  and  Jackson,  1976;  Kaula  et  al.,  1975), 
as  well  as  for  the  geodetic  problem  of  determining  spherical  har- 
monics of  the  geopotential  from  satellite  observations  by  least- 
squares  collocation  (Schwarz,  1974,  1975). 

A general  solution  of  the  underdetermi ned  system  (2-1) 
(assumed  full  rank)  is 

£ = CBt(BCBt)~1x  , (2-2) 

where  £ is  a N x N matrix  such  that  BCBT  is  a regular  n x n 
matrix.  Otherwise  £ is  arbitraty:  different  solutions  are  ob- 
tained by  different  choices  of  £ , and  this  gives  an  infinite 
set  of  possible  solution  vectors  £ . 

It  is  immediately  verified  that  (2-2)  satisfies  the  given 
system  (2-1).  Less  obvious  but  also  well  known  from  the  theory  of 
generalized  inverses  (cf.  Bjerhammar,  1973  , p . 1 1 6 ) is  the  fact 
that  the  solution  (2-2)  satisfies  the  minimum  condition 

£!£  *£  = minimum  , (2-3) 


-l 


provided  the  inverse  matrix  C 


exists  in  an  appropriate  sense 
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(for  N = * this  implies  convergence  of  any  occurring  infinite 
sums)  . 

Usually  the  measurements  x will  be  affected  by  unknown 
observational  errors,  denoted  by  n , the  notation  follows  the 
terminology  of  time  series;  "_s"  stands  for  signal,  and  "_n" 
stands  for  “noise".  (There  is  hardly  any  danger  of  confusing  n , 
the  noise  vector,  with  n , the  number  of  observations.)  Then 
(2-1)  is  to  be  replaced  by 

£s  + n = x . (2-4) 

An  appropriate  minimum  condition,  instead  of  (2-3),  is  now 

s_1£”is_  + nTj)”iji  = minimum,  (2-5) 

where  £ and  £ are  symmetrical  matrices  that  can  be  interpreted 
statistically  as  covariance  matrices:  £ is  the  covariance  matrix 

of  the  signal  s_  , and  D is  the  covariance  matrix  of  the  noise 
jh  , that  is,  of  the  observational  errors. 

The  solution  of  (2-4)  under  the  minimum  condition  is  found 

to  be 


_s  = CBt(BCBt  + 0)  ix  . (2-6) 

This  solution  has  Deen  used  for  geophysical  inverse  problems 
(Wiggins,  1972  , pp. 260-1  ; Burkhard  and  Jackson,  1976  , p . 1 5 1 4 ) . 

It  is  also  the  solution  obtained  by  least-squares  collo- 
cation (Schwarz,  1974).  This  is  clear  from  the  condition  (2-5) 
and  may  also  be  derived  directly  as  follows. 

The  basic  collocation  model  is 

x = £X  + s/+n  (2-7) 
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(Moritz,  19/2,  p.7),  where  s'  is  the  r a ruioi.i  "signal  |jart"  of 
the  observations  x_  and  X denotes  non-random  (systematic) 
parameters.  If  X_  = 0 , the  model  (2-7)  reduces  to 

x = l'  + jn  . (2-8) 

If  the  vector  _s'  is  expressed  as  a linear  combination  6_s  of 
"basic  signals"  (in  the  geodetic  case,  e.g.,  spherical  harmonic 
coefficients),  then  (2-8)  becomes 

x = Bs  + n , (2-9) 

which,  in  fact,  is  the  model  (2-4),  the  "parameters"  s^  now 
being  treated  as  random  variables. 

The  collocation  solution  follows  from  eq.  (2-36)  of 
(Moritz,  1972,  p.15)  for  X = 0 : 


s 


C C” 1 x 

-sx-xx- 


(2-10) 


which  is  essentially  nothing  but  the  fundamental  Wi ener-Ko 1 mogorov 

prediction  formula  (cf.  Liebelt,  1967;  he  calls  it  Gauss-Markov 

theorem).  The  matrix  C is  the  covariance  matrix  of  the  vector 

x , and  C is  the  cross-covariance  matrix  between  the  vectors 
— — sx 

*■§.  and  x.  • 

We  assume  the  vectors  s^  and  ji  to  be  uncorrel  ated . 

Then  the  application  of  covariance  propagation  (Moritz,  1972,  p.97) 
gives  readily. 

Cxx  = BCBT  + 0 , (2-11) 


I 


(2-12) 


£ being  the  covariance  matrix  of  the  vector  s . Inus,  in  our 
case,  (2-10)  in  fact  becomes  (2-6). 

Let  us  now  turn  to  the  accuracy  of  the  estimated  signal 
s_  . This  accuracy  is  usually  defined  by  the  error  covar  i ance 
mat r i x consisting  of  error  variances  (squares  of  standard 

errors)  as  diagonal  terms  and  error  covariances  as  off-diagonal 
terms. 

If  we  abbreviate  (2-6)  as 


s = Lx  , 


(2-13 


then  £ss  is  given  by  eq.  (3-20)  of  (Moritz,  1972,  p.29): 


E = C - LCr 

— ss  — ss  — sx 


C L 1 + LC  L 

S X XX" 


(2-14, 


With  £ ^ = £ and  (2-11)  and  (2-12),  this  becomes 
£ss  = £ - L3£  - £BT£‘+  L(8CBT+0)LT  . 


This  is  readily  given  the  form 


(2-15) 


w n e r e 


L D 
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2 denoting  the  unit  matrix.  This  form  provides  an  elegant  de- 
composition of  the  error  covariance  matrix  which  has  the 

following  interpretation  (Burkhard  and  Jackson,  1976  , p . 1514)  . 
Let  us  write  (2-13)  in  the  form 

£ = Lx  , (2-18) 

where  the  notation  £ is  to  indicate  that  we  are  dealing  wu.i 
the  estimated  value  of  the  signal.  We  now  substitute  into  this 
f ormu la  eq . ( 2 -4 ) , 

x = Bs  + n , (2-19) 

in  which  s^  and  jn  denote  the  "true"  values.  The  result  is 

s = LBs  + £n  , ( 2 - 20 ) 


so  that  the  "true  error"  of  the  signal,  defined  as  estimated 
value  s minus  true  value  s , is  given  by 

£ - s = ( L_B  - £)s_  + L_n  . (2-21) 

The  first  term  on  the  right-hand  side, 

ii  = (LB  - £)s  , (2-22) 

denotes  the  rescl vi ng  error , due  to  the  deviation  of  the  product 
L B from  the  unit  matrix  £ : if  the  system  (2-1)  could  be 
exactly  solved--^  being  a regular  square  matrix--then  e would 
be  zero.  The  second  term, 


Ln  , 


(2-23) 
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expresses  simply  the  effect  of  data  errors  propagating  into  the 
solution.  Clearly,  £ is  the  covariance  matrix  of  e,  and  £ 
the  covariance  matrix  of  e , E and  E simply  add  in  (2-15) 
because  e_.  and  _e_,  are  u ncor  rei  ated , Deing  linear  functions 
of  the  uncorreiated  vectors  s_  and  _n  , respectively. 

We  have  at  length  considered  underdetermi ned  systems,  be- 
cause they  are  typical  for  geophysical  inversion  problems  and 
for  the  determination  of  the  gravitational  field.  As  regards 
(full-rank)  overdeterml ned  systems,  the  model  (2-4)  can  still  be 
used  with  8 as  a "standing"  instead  of  a "lying"  rectangular 
matrix,  the  vector  s_  having  now  less  components  than  the  vector 
x_  , N < n . Tne  condition  (2-5)  gives  again  a solution  which  is 
formally  identical  to  (2-6). 

The  new  feature,  peculiar  to  the  overdetermi ned  case 
n > N , is  that  now  the  following  well-known  transf ormat i on 
(cf.  Liebelt,  1967,  p.30)  can  be  applied: 

CBT(BCBT+D)~1  = (BV1^-1)-1^-1  . (2-24) 

Thus,  for  the  overdetermined  case,  (2-6)  is  equivalent  to 

1 = (iT^'1l+£"1)'1IT£~12i  • (2-25) 

This  form  has  the  advantage  of  permitting  an  easy  transition  to 
least-squares  adjustment  by  parameters:  the  condition  for  this 
case, 


n 1 D 1 n = mi n i mum  , (2-26) 

arises  from  (2-5)  by  letting  ' £~ 1 * 0 ; in  the  same  way,  (2-25) 
becomes 


r 
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s = (BtD  1B)"1BtD  ‘lx  , (2-27 ) 

which  is  the  classical  least-squares  solution  of  (2-4)  if  it  is 
an  overdetenni ned  system. 

It  may  be  mentioned  that  the  two  elementary  cases  of 
classical  least-squares  adjustment,  namely  adjustment  by  para- 
meters and  adjustment  by  conditions,  come  out  as  limiting  cases 
of  the  general  model  (2-6).  They  are,  so  to  speak,  at  extreme 
ends  of  the  scale:  for  overdetermined  systems,  (2-6)  gives  (2-27) 
for  C_1  -*•  0 , as  we  have  just  seen.  For  u nderdetermi  ned  systems, 
(2-6)  becomes  (2-2)  on  letting  D.  -+  0 . However,  (2-6)  is  for- 
mally identical  to  the  result  of  least-squares  adjustment  by 
conditions,  on  interpreting  s^  now  as  measuring  errors;  cf.  eq. 
(12.1:11)  of ( Bj er hammar , 197  3 , p . 1 6 6 ) with  our  eq.  (2-2)  for 
C = l . 

Complete  Sets  of  Solutions.-  An  interesting  question  is 
the  problem  of  determining  the  set  of  a 1 1 “reasonable"  gravi- 
tational fields  that  are  compatible  with  the  set  of  all  avail- 
able geodetic  data.  Each  measurement  (such  as  terrestrial  gravity 
or  satellite  data)  gives  an  equation,  which,  on  linearization, 
provides  one  of  the  linear  equations  of  a system  (2-1);  we  assume 
n independent  observations.  The  unknowns  s^  are  to  be  the  set 
of  all  zonal  and  tesseral  harmonics  in  the  anomalous  gravi- 
tational field;  thus  N = » . 

For  simplicity  we  assume  errorless  data  (n_  = 0)  ; thus 
the  equations  (2-1)  have  to  be  satisfied  exactly.  As  we  have 
mentioned,  such  a complete  set  of  solutions  is  obtained  by  let- 
ting the  matrix  C vary  in  the  solution  (2-2).  However,  this 
representation  is  computationally  inconvenient  since  each  time 
a large  matrix  BC B1  would  have  to  be  inverted  anew. 

Practically  more  suitable  is  the  representa t i on , well- 
known  from  the  theory  of  singular  matrix  inverses  (cf.  Bjerhammar, 
1973  , p . 378)  , 
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H”x  + ( 1 H t<)u  , ( 2 28) 

w h e r e 

B'  = CoBt(BC  Bt)"1  (2-29) 

is  a fixed  generalized  inverse  of  j5  , to  be  computed  once  for  all 
with  an  assumed  fixed  matrix  C_c  . The  vector  u_  is  arbitrary; 
the  desired  set  of  solutions  is  obtained  by  letting  u_  run  through 
the  set  of  ail  possible  vectors  of  infinitely  many  components  for 
which  the  infinite  sums  converge. 

Having  thus  obtained  the  complete  set  of  solutions  (2-28), 
one  may  select  special  solutions  by  imposing  suitable  conditions. 

The  geodetically  most  meaningful  condition  is  probably  optimal 
accuracy  of  the  result,  expressed  in  terms  of  minimum  error  vari- 
ances. This  condition  is  known  to  be  equivalent,  for  errorless  data, 
to  the  condition  (2-3),  C being  the  signal  covariance  matrix 
(cf.  Moritz,  1972,  pp.38  and  122).  Under  this  condition  the  solution 
is  given  by  (2-2),  For  this  case  we  may  put 

u = CBt(BCBt)"1  - B"  x ; (2-30) 

in  fact, on  substituting  this  expression,  (2-28)  becomes  (2-2). 

A geophys i ca 1 1 y interesting  side  condition  is  that  the 
solution  ^ belong  to  a convex  set  described  by  inequalities  such 
a s 

s > 0 , 2s,  + s = 5 , etc. 

Problems  of  minimizing  the  error  variance  (or  other  functions) 
subject  to  such  side  conditions  are  dealt  with  in  optimization 


L 
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theory  (linear  and  nonlinear  prog ramm i ng ) . 

In  the  problem  of  determining  mass  structures  from  gravity 
data,  the  requirement  that  the  obtained  density  be  nonnegative 
would  lead  to  such  side  conditions. 

In  the  geodetic  inversion  problem,  the  determination  of  i 

the  external  gravitational  field,  such  conditions  have  never  been 
used  and  do  not  seem,  in  general,  to  have  practical  importance. 

However,  the  condition  that  the  total  gravitational  field  should 
be  generated  by  masses  that  are  everywhere  nonnegative,  is  a 
physical  requirement  to  be  met.  In  "reasonable"  solutions  it 
seems  to  be  more  or  less  automatically  satisfied,  but  a closer 
look  on  the  convexity  problem  in  geodetical  applications  may  not 
be  without  theoretical  interest,  as  has  been  pointed  out  to  the 
author  recently  by  Prof.  P.C.  Sabatier  and  Prof.  0 . D . Jackson. 

We  finally  mention  the  beautiful  geometrical  treatment  of 
underdetermined  and  overdetermined  systems  of  linear  equations  in 
(Lanczos,  1964,  chapter  3). 

3 . Functional  Representation 
of  the  Gravitational  Field 

Let  us  now  turn  to  an  apparently  quite  different  problem, 
the  representation  of  the  earth's  external  gravitational  field  by 
appropriate  analytical  functions. 

Let  the  anomalous  potential  T be  approximated  by  a linear 
combination  f of  suitable  base  functions  <t>  i *4  2 3 » • • • : 

T(P)  = f(P)  = l b <j>  ( P ) , (3-1) 


P denoting  the  space  point  at  which  these  functions  are  being 
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considered,  and  denoting  suitable  coefficients. 

Since  T is  harmonic  outside  the  earth's  surface,  all  base 
functions  must  be  harmonic  functions,  too.  They  may,  for 

instance,  be  spherical  harmonics  or  the  potentials  of  point  masses 
suitably  distributed  below  the  earth's  surface. 

How  are  the  coefficients  b^  to  be  determined?  Assume, 
for  the  sake  of  simplicity  and  definiteness,  that  we  are  given 
errorless  values  of  T at  n space  points  P..  . Then  it  is 
reasonable  to  use  the  approximation  of  T by  a linear  combination 
(3-1)  of  n base  functions  and  to  postulate  that  f exactly 
reproduces  T at  the  n given  points.  Putting 


TCP,)  = f(Pi)  = f±  , 

we  thus  have  the  conditions 


JAW  ■ fi 


(3-2) 


(3-3) 


which  are  n linear  equations  for  the  n unknowns  bR  , which 
can  be  solved  provided  they  are  linearly  independent. 

With 


Mpi) 


AiX 


we  thus  have 


y a 


k = 1 


b 

i k k 


f . 


or  in  matrix  notation, 


(3-4) 


(3-5) 


Ab 


jC 


(3-6) 


14 


with  the  solution 


b = A~ 1 f_  . (3-7) 

Thus  the  determination  of  the  coefficients  always  involves 

the  solution  of  a n x n system  of  linear  equations,  or  the  in- 
version of  a n x n matrix. 

An  exception  would  be  the  case  that  the  matrix  A re- 
duces to  the  unit  matrix,  that  is 


1 , i = k 

i 

) 

,0  , i i k 


(3-8) 


the  base  functions  for  this  case  being  denoted  by  x(P)  instead 
of  <j>(P)  . This  means  that  the  functions  xk  are  zero  at  all  data 
points  except  one,  at  which  they  assume  the  value  1 . Such 
functions  are  called  sample  functions;  they  have  the  pleasant 
property  that,  by  (3-7),  b = f , so  that  (3-1)  reduces  to 


f(P)  = l fkxk(P)  . (3_9) 

*=1 

\ 

the  coefficients  being  simply  the  given  values  of  the  function 
at  the  data  points  PR  . 

Of  course,  it  cannot  be  expected  in  general  that  A in 
fact  reduces  to  the  unit  matrix,  but  even  if  this  is  not  the  case, 
we  can  associate  a set  of  sample  functions  xk  to  the  given  base 
functions  , by  putting 

xk(P)  - 


(3-10) 
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where  A(  l)  denote  the  elements  of  the  matrix  A- 1 inverse 

J K — 

to  A . In  fact,  putting  P = P.  in  (3-10)  we  get 


Xk(P1)=.IiA‘-1)Ai. 


: l a.  ,a(::) 

j=l  13  3k 


= 5 


ik  ’ 


(3-111 


according  to  the  definition  of  the  inverse  matrix,  so  that  (3-8) 
is  satisfied.  Furthermore,  the  substitution  of  (3-10)  into  (3-9) 
gives 


f(P)  = l i Aj<k1)vj(p) 

j-lk-i  3 3 


n 

= l b ♦ (P)  (3-12) 

3=1  D D 

by  (3-7),  so  that  the  sample  function  development  (3-9)  is  iden- 
tical to  the  original  representat i on  (3-1). 

In  this  way  we  can,  indeed,  associate  to  each  system  of 
base  functions  ....  and  to  each  configuration  of  points 

P1,P2 Pj(  a system  of  sample  functions  x1»x2»-**»x  • This 

has  advantages  if  the  same  point  configuration  is  used  with 
different  data  sets;  then  only  the  f^  change,  but  the  sample 
functions  xk(P)  remain  the  same.  As  a matter  of  fact,  also  this 
approach  requires,  according  to  (3-10),  the  inversion  of  the 
n x n matrix  A , but  it  need  to  be  performed  only  once  for  a 
given  point  configuration. 

The  type  of  sample  functions  best  known  in  geodesy  are 
the  functions  of  Giacaglia  and  Lundquist  (1972)  based  on  the 
system  of  spherical  harmonics  for  regular  point  co nf i gu ra t i ons 
on  the  sphere.  For  sucn  point  configurations 


J 
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it  presents  computational  advantages  and  is  precisely  equivalent 
to  a spherical -harmonic  representation  up  to  a certain  degree. 

Spherical  harmonics  and  sample  functions  based  on  them 
are  suitable  for  a global  representat i on . For  local  represen- 
tations, potentials  of  buried  masses  or  related  functions  (Dufour 
and  Kovalevsky,  1970)  are  more  suitable;  another  possibility  is 
least-squares  interpolation  to  be  considered  now. 

The  approximation  of  the  potential  T by  a finite  llr.ear 
combination  of  n base  functions  $k  can,  of  course,  reproduce 
T only  at  n points,  the  data  points.  At  other  points,  f will 
deviate  from  T , that  is,  the  error 

£p  = T ( P ) - f (P)  (3-13) 

will,  in  general,  differ  from  zero.  Of  considerable  importance  is 
obviously  the  question  whether  there  exists  a system  of  base 
functions  ^(P)  for  which,  at  every  point  P , the  mean  square 
interpolation  error  mp  , defined  by 

"«p  = M(e2)  (3-14) 

attains  a minimum,  M denoting  a suitably  defined  statistical 
expectation  value. 

The  answer  to  this  question  is  yes  (which  is  by  no  means 
a matter  of  fact).  This  leads  to  least-squares  interpolation  due 
to  N.  Wiener  and  A.N.  Kolmogorov.  A derivation  is  found  in 
(Heiskanen  and  Moritz,  1967,  p.268).  The  result  is 

(p)  = C(P,Pk)  , (3-15) 

where  C ( P , P ) is  tne  covariance  function  of  T between  points 
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P and  Py  . Now  tnis  covariance  function  admits  the  r epr esenta t i on 
(Mor i tz  , 1 9 7 2,  p 88 ) 

C(P,Q)  - )'  kjr ~)  ’’  1 P;j  ( c o s :. ) , (3-16) 

a=o  ’ P Q 

where  k and  R are  constants,  r and  r are  the  spat.  a. 

i,  P Q 

radius  vectors  of  P and  Q , and  P are  Legendre's  polynomials 
as  functions  of  the  angle  p between  rp  and  r^  . This  re- 
presentation shows  that  C ( P , P ) , considered  as  a function  of  P , 
is  a harmonic  function,  analytic  everywhere  outside  a certain 
sphere.  If  the  radius  of  this  sphere  is  chosen  so  that  the  sphere 
lies  completely  inside  the  earth,  then  the  functions  (3-15)  will 
be  analytic  everywhere  outside  and  on  the  earth's  surface  and  can, 
therefore,  oe  used  for  representing  the  external  gravitational 
potential. 

It  is  a general  principle  of  least. - squares  estimation  that 
the  result  does  not  depend  strongly  on  the  a priori  covariances  i 

used.  Therefore,  if  the  actual  covariance  function  (3-16)  is  a 
function  which  is  too  complicated  for  practical  use,  one  may, 
witnout  appreciable  loss  of  accuracy,  replace  the  actual  function 
by  a suitable  approximation  of  a simpler  and  more  manageable 
analytical  form  (in  practice,  one  will  anyway  fit  a rather  simple 
analytical  expression  to  empirically  obtained  covariance  data). 

’’herefore,  also  the  functions  (3-15)  will  be  precisely  defined 
harmonic  functions  of  a relatively  simple  analytical  form. 

By  (3-4)  and  (3-15)  we  have  now 

A ik  = c(Pi,pk)d=fC.k  . (3-17) 

Defining  the  matrix  i by 


C = (C . J = A , (3-18; 

— i k — v 


T 
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we  get  from  (.<-/) 


(3-19) 


On  introducing  the  vector 

C(P)  = : C ( P , P , ) , C(P,P2)  ,.  . . , C(P,Pn)  , 


(3-20) 


(3-1)  thus  takes  the  form 

f(P)  = C(P)C_1f  , (3-21) 

which  is  again  the  Wiener-Kolmogorov  prediction  formula  (2-10). 

It  is  clear  that,  on  the  basis  of  the  functions  (3-15) 
for  least-squares  prediction,  corresponding  sample  functions  can 
be  defined  by  (3-10),  each  of  which  is  zero  at  all  data  points 
except  one,  at  which  it  takes  the  value  1. 

Least-squares  prediction  shares  with  all  representation 
methods  of  type  (3-1)  the  disadvantage  that  a (usually  large) 
n x n matrix  has  to  be  inverted;  it  has,  however,  the  advantage 
that  in  this  particular  case  the  matrix  to  be  inverted  is  symmetr i c , 
being  a covariance  matrix  (in  the  general  case,  the  matrix  A 
is  not  symmetr i c ) . 

Taken  in  the  form  just  described,  the  application  of  least- 
squares  prediction  in  geodesy  is  rather  limited,  essentially  to 
the  prediction  of  gravity  anomalies.  In  fact,  however,  geodetic 
data  are  of  many  different  types:  horizontal  and  vertical  angles, 
distances,  gravity  measurements,  deflections  of  the  vertical, 
satellite  observations  of  many  kinds,  gravity  gradients,  etc.  All 
these  geodetic  data  obviously  share  the  properties  that  they 
depend  on  the  earth's  gravitational  field;  therefore,  on  lineari- 
zation, all  of  them  may  be  expressed  a linear  functionals  of  the 
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anomalous  potential  T . (Linear  Functionals  are  by  no  means 
restricted  to  linear  functions  but  they  include  differential  and 
integral  formulas  such  as  Stokes'  integral.) 

So  more  general  and  geodetically  more  important  is 
the  problem  to  fit  the  representation  (3-1)  to  the  data,  so  that 
the  n given  functionals  of  T are  exactly  reproduced.  This  is 
the  principle  of  collocation,  which  is  frequently  used  in  numer- 
ical mathematics  (cf.  Collatz,  1966).  If  the  base  functions 
? ( P ) are  again  to  be  determined  by  the  least-squares  condition 
of  minimum  interpolation  error  at  all  points  P , then  it  is 
natural  to  call  the  method  least-squares  collocation. 

The  solution  has  again  the  same  formal  structure  as  in 
least-squares  prediction.  The  result  is 

?k(P)  = C(P,xk)  , (3-22) 

which  is  the  covariance  between  T(P)  and  the  measurement  ; 

it  is  a function  of  the  point  P that  is  again  harmonic  and 
analytic.  The  coefficients  bk  forming  the  vector  b^  are  deter- 
mined by 

- = -xx-  ’ (3-23) 

where  C is  the  autocovariance  matrix  of  the  observation 

-XX 

vector  x = (x  ) . 

— k. 

Thus,  with  errorless  data,  least-squares  collocation 
determines  the  analytical  form  of  the  functions  <f>k  by  the  re- 
quirement of  optimal  accuracy,  whereas  the  data  are  exactly  repro- 
duced. 

In  general,  the  geodetic  data  will  be  affected  by  random  , 

measuring  errors.  Also  in  this  case  the  condition  m = minimum 
can  be  applied  and  determines  simultaneously 

A 
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(1)  the  best  analytical  expression  for  the  functions 
and 

(2)  the  best  values  for  the  coefficients  b 

K 

Even  in  this  case,  the  formal  expressions  are  the  same 
as  before,  again  being  given  by  (3-22)  and  (3-23).  The  measuring 
errors  have  no  influence  on  the  choice  of  4>  k by  (3-22),  so 
that  ; again  represent  pure  gravitational  covariance  functions, 
that  is,  analytical  and  harmonic  functions;  again  the  least-squares 
principle  serves  only  to  single  out  the  most  suitable  analytical 
expression  for  the  base  functions  $.  among  the  many  possible 
choices. 

Where  statistics  enters  directly  is  the  determination  of 
the  coefficients  b , which  is  done  in  such  a way  that  the  effect 

K 

of  random  data  errors  is  minimized  (therefore,  the  data  are  no 
longer  exactly  reproduced);  in  statistical  terminology,  we  nave 
a "best  linear  estimate": an  unbiased  estimate  of  minimum  vari- 
ance. 

Expressions  analogous  to  (3-1)  may  be  given  for  any  other 
quantity  of  the  anomalous  gravitational  field  (called  "signal"), 
such  as  geoidal  heights,  deflections  of  the  vertical,  gravity 
anomalies,  etc.  The  coefficients  b^  remain  the  same  since  they 
depend  only  on  the  data  x by  (3-23).  What  changes  are  the  base 
functions  c k ; the  new  base  functions  are  derived  by  simple 
analytical  operations  such  as  differentiation,  since  a linear 
operation  performed  on  (3-1)  acts  on  the  base  functions  <{)..  only. 

Since  these  base  functions  are  covariance  functions,  these 
operations  are  special  cases  of  covariance  propagation  which  plays 
a basic  role  in  least-squares  collocation:  it  has  to  carry  the 
precise  mathematical  structure  of  the  gravitational  field  (cf. 
Moritz,  1972,  pp. 94-99). 

Inclusion  of  Systematic  Parameters .-  A last  restriction 
has  to  be  removed  before  least-squares  collocation  can  be  fully 
applied  to  general  geodetic  problems.  So  far  we  have  assumed  that 


we  deal  only  with  quantities  that  have  zero  statistical  expectation 
such  as  the  elements  of  tne  anomalous  gravitational  field,  so 
that  systematic  effects  had  to  be  removed  beforehand;  we  shall 
now  free  ourselves  from  this  restriction. 

We  shall,  therefore,  consider  the  following  model: 

x = AX  + s + n , (3-24) 

where  the  vector  .x  comprises  the  measured  quantities,  s^ 
representing  the  part  due  the  anomalous  gravitational  field  and 
jn  denoting  random  measuring  errors.  The  new  component  is  A_X  , 
where  the  vector  X^  comprises  systematic,  nonrandom  parameters 
and  A is  a given  matrix  of  coefficients. 

This  model  is  general  enough  to  encompass  all  conceivable 
geodetic  measurements.  In  fact,  any  geodetic  measurement  may  be 
split  up,  according  to  (3-24),  into  three  parts: 

1.  A systematic  part  AX^  involving,  on  the  one  hand,  the 
ellipsoidal  reference  system  and,  on  the  other  hand,  other  para- 
meters and  systematic  errors  (originally  non-linear  functions  are 
again  thought  to  have  been  linearized  by  Taylor's  theorem); 

2.  A random  part  s^  (of  zero  expectation)  expressing  the 
effect  of  the  anomalous  gravity  field,  and 

3.  Random  measuring  errors  n,  . 

As  an  example,  consider  a measurement  of  gravity,  g . 

Here  AX^  represents  normal  gravity  y , as  well  as  systematic 
errors  such  as  gravimeter  drift;  js  is  the  gravity  anomaly  pg  ; 
and  ji  stands  for  the  measuring  error.  Other  examples  will  be 
found  in  (Moritz,  1972,  pp. 70-76). 

The  formulas  for  estimating  X^  and  s^  may  be  derived 
from  two  different  but  equivalent  minimum  principles: 

1.  From  a least-squares  principle  corresponding  to  (2-5); 

2.  From  the  condition  of  minimum  variance  (least  standard 
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errors  for  estimated  X and  s^  ). 

The  result  is: 

X = ( ATC~ 1 A ) ” 1 ArC - 1 x , 

(3-25) 

s * C C-1(x  - AX)  . 

— — sx  xx  -  4 * *  7 

(3-26) 

The  first  equation  is  analogous  to  classical  1 ea s t - squa re s adjust- 
ment by  parameters,  except  that  the  covariance  matrix  £xx  in- 
cludes now  covariances  of  the  signal  as  well  as  those  of  the 
measuring  errors.  The  second  equation  is  a fairly  obvious  gener- 
alization of  the  Kol mogorov -Wi ener  formula  (2-10)  or  (3-21)  to 
the  case  in  which  the  expectation  of  x is  AX  rather  than  zero. 

These  formulas  are  an  extension  of  tlie  corresponding  re- 
sults for  time  series  (Grenander  and  Rosenblatt,  1957,  p.87). 

The  present  method,  least-squares  collocation  with  para- 
meters, may  be  considered  as  a combination  of  least  squares  adjust- 
ment and  least-squares  prediction  into  a unified  scheme,  which 
makes  possible  the  use  of  all  conceivable  geodetic  da ta --c 1 a s s i ca 1 
angle  and  distance  measurements,  gravity  measurements,  satellite 
data  of  different  kind,  etc. --to  obtain  the  geometric  position  of 
points  of  the  earth's  surface  as  well  as  the  gravitational  field. 

4 . The  Many  Facets  of  Collocation 

The  Nature  of  Least-Squares  Collocation.-  In  sec.  2 we 

have  seen  that,  in  important  cases,  least-squares  collocation 

reduces  to  a least-squares  solution  of  an  underdetermined  system 
of  linear  equations.  It  might  thus  seem  that  1 ea s t - squa re s collo- 
cation is  essentially  nothing  else  than  classical  least-squares 
methods  known  from  adjustment  computations. 


In  trying  to  answer  this  question,  we  shall  first  consider 
"simple"  least-squares  collocation,  that  is  without  estimation 
of  systematic  parameters,  so  that  we  have  A = 0 in  (3-24)  and 
(3-26);  this  was  also  done  in  sec. 2.  Physically  this  means  that 
the  observations  _x  are  all  quantities  of  the  anomalous  gravity 
field. 

In  fact,  every  problem  of  simple  least-squares  collocation 
can,  in  principle,  be  formulated  as  a system  of  linear  equations 
of  the  form  discussed  in  sec. 2.  Every  quantity  of  the  anomalous 
gravity  field,  for  instance  a geoidal  height,  a gravity  anomaly, 
or  a deflection  of  the  vertical,  can  be  expressed  as  an  infinite 
series  of  spherical  harmonics,  the  coefficients  of  which  constitute 
the  unknowns,  the  number  N of  which  is  obviously  infinite.  (Con- 
vergence problems  may  be  overcome  by  the  use  of  Runge's  theorem 
(Krarup,  1969,  p.54;  Moritz,  1971,  p.79).)  Thus  we  always  get  a 
system  of  equations  of  type  (2-9),  as  outlined  in  sec. 2. 

If  N were  finite,  then  the  solution  (2-6)  would,  in  fact, 
correspond  to  classical  1 ea st -squares  since  the  condition  (2-5)  can 
be  written  in  the  form 

vTPv  =minimum  (4-1) 

with 

| C'1  0 i 

1 

P = | 

! 0 D"1  (4-2 ) 

L j 

and  (2-4)  takes  the  form  of  a condition  equation  for  the  vector  v . 

What  difference  does  the  fact  N = * make?  The  signal  space 
is  now  no  longer  finite-dimensional  but  it  is  infinitely-dimensional 
Hilbert  space.  The  elements  of  Hilbert  spaces  may  be  infinite 


vectors,  such  as  s_  , but  also  functions  of  a certain  kind.  In 
the  gravitational  case,  we  conveniently  consider  spaces  of  har- 
monic functions,  more  precisely,  of  functions  harmonic  outside 
a certain  sphere. 

There  is  a one-to-one  correspondence  between  such  func- 
tions and  infinite  vectors  _s  (of  norm  s_1' _s  < ) by  taking  £ 

as  the  vector  of  coefficients  of  the  spher i ca 1 -ha rmon i c expansion. 
This  makes  it  possible  to  avoid  the  use  of  infinite  vectors  (ana 
correspond  i ng  convergence  problems)  by  working  with  functicr... 
and  linear  operations  on  them.  In  this  way,  formulas  of  Kolmogorov- 
Wiener  type  (3-21)  are  obtained,  which  operate  with  finite 
n x n matrices. 

Another  way  of  linking  collocation  (including  systematic 
parameters)  with  adjustment,  this  time  in  f i n i te -d i mens i ona 1 
space,  was  given  in  (Moritz,  1972,  pp.  12-15).  There,  however,  the 
signal  values  to  be  computed  do  not  enter  in  the  condition  equa- 
tions although  they  do  enter  into  the  minimum  principle  y_J  Pv_  . 

Thus  the  condition  equations  do  not  provide  a complete  formulation 
as  they  do  in  adjustment  computations,  but  must  be  supplemented 
by  additional  considerations:  the  new  signals  are  related  to  the 
observations  not  by  the  condition  equations,  but  through  joint 
covariances.  (The  importance  of  joint  covariances  is  well  known 
from  wide-sense  stationary  stochastic  processes  where  they  carry 
the  total  statistical  structure.) 

So  the  relation  of  collocation  and  adjustment  is  very 
elusive:  the  analogies  are  striking,  but  at  the  very  moment  when 
we  think  that  we  have  hit  on  an  exact  identity  we  must  recognize 
a difference  in  a fine  but  essential  point.  For  a more  detailed 
discussion  see  also  (Rummel,  1976). 
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A much  closer  relationship  is  between  least-squares  collo- 
cation and  the  theory  of  prediction  of  stationary  stochastic  pro- 
cesses. In  fact,  the  anomalous  gravity  potential  may  be  considered 
as  a stochastic  process  on  the  sphere,  or  rather  as  a spatial 
stochastic  process  that  is  harmonic  outside  a sphere. 

The  theory  of  stochastic  processes  provides  a very  con- 
venient mathematical  formalism  and  a statistical  i nter pr eta t i on 
and  terminology  (e.g.,  covariance  functions).  In  fact,  we  have 
seen  at  the  end  of  the  preceding  section  that  the  basic  collocation 
formulas  (3-25)  and  (3-26)  are  precise  analogues  of  the  corre- 
sponding results  for  prediction  and  parameter  estimation  of  time 
series. 

Some  kind  of  statistical  interprediction  is  desirable  to 
get  the  geodetically  important  concept  of  accuracy  (standard  errors 
of  prediction)  into  the  picture.  How  seriously  the  statistical 
i nterpretat i on  is  taken,  is  a matter  of  controversy  and  also  of 
personal  taste.  The  present  author  (Moritz,  1972,  secs. 8 and  9) 
favors  an  interpretation  in  terms  of  Norbert  Wieners  "covariance 
of  individual  functions"  to  take  into  account  the  fact  that  there 

is  only  one  earth  and  to  avoid  difficulties  associated  with  the 

stochastic  process  interpretation  as  pointed  out  by  Lauritzen 
(1973):  it  is  impossible  to  find  a stochastic  process,  harmonic 

outside  the  sphere,  that  is  both  Gaussian  and  ergodic. 

It  is  possible  largely  to  play  down  tne  statistical  aspects, 
emphasizing  Hilbert  space  geometry  and  considering  the  covariance 
function  as  a kernel  function  in  Hilbert  space,  as  Krarup  did  in 
his  fundamental  paper  (1969).  In  fact,  an  interpolation  formula 
formally  identical  to  the  Kolmogorov-Wiener  prediction  formula  is 
obtained  in  the  theory  of  kernel  functions  in  Hilbert  space 
(Meschkowski , 1962,  p . 114)  in  a completely  "deterministic"  way, 
without  using  any  statistics. 

It  is  well  to  emphasize  the  analytical  nature  of  kernel 
f u nc t i ons - -wh i c h are  precisely  defined  narmonic  f unct i ons--to 
avoid  the  wrong  impression  that  1 east-squares  collocation  "messes 
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up  everything  statistically"  and  to  show  that  we  have  an  under- 
lying completely  "clean"  analytical  model.  The  elementary  rea- 
soning in  sec. 3 of  the  present  paper  is  intended  to  serve  the 
same  purpose. 

Why,  then,  use  a special  name,  "least-squares  collocation" 
for  the  present  geodetic  method  and  not  simply  call  it  prediction? 
The  term  prediction  is  usually  understood  to  comprise  interpo- 
lation and  extrapolation  of  time  series  and  data  of  the  same 
kind,  for  instance,  gravity  anomalies.  The  essence  of  the  present 
method,  however,  is  the  use  of  heterogeneous  data,  which  are  linear 
functional s of  the  anomalous  gravitational  potential.  As  we  have 
seen  in  the  preceding  section,  the  fitting  of  a function  to  given 
functionals  is  precisely  the  feature  of  collocation  as  understood 
in  numerical  mathematics. 

The  name  "least-squares  collocation"  thus  seems  to  be  quite 
suitable  to  characterize  the  present  method  as  a representation 
of  the  anomalous  gravitational  field  by  "clean"  analytical  func- 
tions, which  are  selected  and  the  coefficients  of  which  are  deter- 
mined by  a least-squares  principle. 


■n  ■ 
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